library(R6)
suppressMessages(library(tidyverse)); suppressMessages(library(magrittr))

# v1
ml <- R6::R6Class(
  classname = "ml",
  public = list(
    exprs = NULL,
    pheno = NULL,
    od = NULL,
    initialize = \(exprs = NULL, pheno = NULL, od = NULL){
      self$exprs <- exprs
      self$pheno <- pheno
      self$od <- od

      index <- list(self$exprs, self$pheno, self$od)

      if(any(map_lgl(index,is.null))){
        features <- c('expr',"pheno","od")[which(map_lgl(index,is.null))]
        cli::cli_alert_warning("缺少{features}信息，需要在`new()`中重新定义")
      }
    },
    help = \(...){
      library(cli)
      ul <- cli_ul()
      cli::cli_text("{.strong 机器学习包装测试:}")
      cli_li("构建新对象需要两个变量，一个是包含选用特征的类似表达谱文件，另一个是分组信息文件。")
      cli_li("随机森林中，包含 w, h, ntree, var_name, top_feature, col_vector等其他参数；")
      cli_li("lasso中包含seed参数;")
      cli_li("svm-rfe默认使用五倍交叉验证；")
      cli_li("使用示例：my_ml_obj$run('lasso')")
      cli_end(ul)

      # message("
      # 机器学习包装测试
      # - 随机森林中，包含 w, h, ntree, var_name, top_feature, col_vector等其他参数；
      # - lasso中包含seed参数;
      # - svm-rfe默认使用五倍交叉验证；
      # - 使用示例：my_ml_obj$run('lasso')
      # ")
    },
    run = function(method = NULL, ...) {
      if(is.null(method)){
        cli::cli_alert_warning("请选择正确的参数，程序暂时不运行任何分析，参考{.arg method = 'lasso|rf|svm_rfe}'")
        help()
        return()
      }

      y <- match.arg(method, c("lasso", "rf", "svm_rfe"))

      f <- switch(y,
        lasso = \(exprs = self$exprs, od = self$od, pheno = self$pheno, seed = 1110, ...){
          # lasso 建模
          library(glmnet)
          set.seed(seed)
          if (any(is.na(exprs[[2]]))) {
            message(cli::style_bold(">> 表达谱有NA，自动去除..."))
            exprs <- na.omit(exprs)
          }

          x_tmp <- exprs %>%
            t() %>%
            as.data.frame() %>%
            rownames_to_column("sample")
          dat <- pheno %>% inner_join(x_tmp)

          x <- dat %>%
            dplyr::select(-any_of(colnames(pheno))) %>%
            as.matrix()
          y <- dat %>%
            dplyr::select(any_of(colnames(pheno))) %>%
            select(-sample) %>%
            pull(1)
          fit1 <<- cv.glmnet(x, y, family = "binomial", type.measure = "deviance", alpha = 1)
          # 绘制图片1
          dir.create(file.path(od, "lasso/"), recursive = T, showWarnings = F)
          pdf(file = file.path(od, "lasso/lasso_fit.pdf"), width = 5, height = 4)
          plot(fit1)
          dev.off()
          # 查看模型各变量系数
          fit1_coef_lambda.min <- coef(fit1, s = fit1$lambda.min)
          # 查看系数非0变量
          fit1.min <- fit1_coef_lambda.min[which(fit1_coef_lambda.min != 0), ]
          # 矩阵转换，使上面结果更方便查看
          fit1.min_2 <- matrix(fit1.min, length(fit1.min), 1)
          rownames(fit1.min_2) <- names(fit1.min)
          colnames(fit1.min_2) <- c("coef")
          lasso_res <- fit1.min_2 %>%
            as.data.frame() %>%
            rownames_to_column("gene")
          lasso_res <- lasso_res[-1, ]
          # (Intercept) -98.7319020

          # 绘制图片2
          fit2 <<- glmnet(x, y, family = "binomial")

          pdf(file = file.path(od, "lasso/lasso_cv.fit.pdf"), width = 5, height = 4)
          plot(fit2, xvar = "lambda", label = TRUE)
          dev.off()

          lasso_fits <- list(x = x, y = y, fit = fit2, cv.fit = fit1)
          saveRDS(lasso_fits, file.path(od, "lasso/lasso_fits.rds"))

          # 输出结果
          write_tsv(lasso_res, file.path(od, "lasso/lasso_res.txt"))

          message(cli::style_bold(">> LASSO binomial 回归分析结束"))
          lasso_fits
        },
        rf = function(
            data = self$exprs, group = self$pheno, od = self$od, w = 8, h = 6.5, ntree = 1000,
            var_name = "random_forest",
            top_feature = 20, col_vector = RColorBrewer::brewer.pal(8, "Set2"), ...) {
          data <- na.omit(data)

          data <- data %>%
            t() %>%
            as.data.frame() %>%
            rownames_to_column("sample") %>%
            inner_join(group)

          group_col <- colnames(group)[2]

          data[[group_col]] <- factor(data[[group_col]])

          dir.create(od, recursive = T, showWarnings = F)

          colnames(data) <- colnames(data) %>%
            str_replace_all("-", "_")
          library(randomForest)
          set.seed(1110)

          mtry_range <- 1:(ncol(data) / 2)

          err_df <- map_df(mtry_range, \(x) {
            set.seed(71)
            rf <- randomForest(
              y = data[[group_col]], x = data %>%
                select(-any_of(group_col)), importance = TRUE, proximity = TRUE, ntree = ntree,
              mtry = x
            )

            err <- mean(rf$err.rate)

            data.frame(mtry = x, mean.error = err)
          })

          mtry_error_p <- err_df %>%
            mutate(type = ifelse(mean.error == min(mean.error), "a", "b")) %>%
            ggplot(aes(x = mtry, y = mean.error)) +
            geom_point(
              shape = 21, aes(color = type),
              size = 2, show.legend = F
            ) +
            geom_line() +
            theme_bw(12) +
            theme(
              legend.title = element_blank(),
              legend.position = "top"
            ) +
            scale_color_manual(values = c("red2", "black")) +
            labs(y = "Mean Error")

          # which.min(err_df[[2]]) 14

          rf <- randomForest(
            y = data[[group_col]], x = data %>%
              select(-any_of(group_col)), importance = TRUE, proximity = TRUE, ntree = ntree,
            mtry = which.min(err_df[[2]])
          )

          error_p <- rf$err.rate %>%
            as.data.frame() %>%
            mutate(trees = seq_len(nrow(.))) %>%
            pivot_longer(-trees, values_to = "Error", names_to = "type") %>%
            ggplot() +
            geom_line(aes(x = trees, y = Error, color = type), size = 1) +
            theme_light() +
            theme(legend.title = element_blank(), legend.position = "top") +
            scale_color_manual(values = col_vector)

          MeanDecreaseAccuracy_p <- importance(rf) %>%
            as.data.frame() %>%
            rownames_to_column("Gene") %>%
            mutate(Gene = str_replace(Gene, "_", "-")) %>%
            arrange(desc(MeanDecreaseAccuracy)) %>%
            .[1:top_feature, ] %>%
            mutate(Gene = factor(Gene, levels = rev(unique(Gene)))) %>%
            ggplot() +
            geom_point(aes(x = MeanDecreaseAccuracy, y = Gene), shape = 1) +
            theme_light() +
            theme(legend.title = element_blank(), legend.position = "top") +
            scale_color_manual(values = col_vector) +
            labs(y = NULL)

          d <- importance(rf) %>%
            as.data.frame() %>%
            rownames_to_column("Gene") %>%
            mutate(Gene = str_replace(Gene, "_", "-")) %>%
            arrange(desc(MeanDecreaseGini))

          MeanDecreaseGini_p <- d %>%
            mutate(Gene = factor(Gene, levels = rev(unique(Gene)))) %>%
            .[1:top_feature, ] %>%
            ggplot() +
            geom_point(aes(x = MeanDecreaseGini, y = Gene), shape = 1) +
            theme_light() +
            theme(legend.title = element_blank(), legend.position = "top") +
            scale_color_manual(values = col_vector) +
            labs(y = NULL)

          saveRDS(rf, file = paste0(od, "/", var_name, "_rf.rds"))

          library(cowplot)

          p1 <- plot_grid(mtry_error_p, error_p, labels = "AUTO", nrow = 2)
          p2 <- plot_grid(MeanDecreaseAccuracy_p, MeanDecreaseGini_p,
            labels = LETTERS[3:4],
            nrow = 1
          )
          p <- plot_grid(p1, p2, rel_widths = c(1, 2), align = "hv")

          ggsave(plot = p, filename = paste0(od, "/", var_name, ".pdf"), width = w, height = h)
          pl <- list(
            mtry_error_p = mtry_error_p, error_p = error_p,
            MeanDecreaseAccuracy_p = MeanDecreaseAccuracy_p, MeanDecreaseGini_p = MeanDecreaseGini_p
          )
          saveRDS(pl,paste0(od, "/", var_name, "_plot_list.rds"))

          return(list('d' = d,'rf_obj' = rf))
        },
        svm_rfe = \(exprs = self$exprs, od = self$od, pheno = self$pheno, ...){
          # SVM
          set.seed(123)
          features <- exprs %>%
            t() %>%
            as_tibble(rownames = "sample")

          labels <- pheno %>% as_tibble()

          exp_t <- inner_join(features, labels, by = "sample") %>%
            column_to_rownames("sample")

          exp_t[is.na(exp_t)] <- 0

          set.seed(1110)
          x <- as.matrix(exp_t[, 1:(ncol(exp_t) - 1)])
          # x <- exp_t
          y <- factor(exp_t[, ncol(exp_t)])

          library(caret)
          control <- rfeControl(functions = caretFuncs, method = "cv", number = 5)
          results <<- rfe(x,
            y = y,
            sizes = seq(0, ncol(x)),
            rfeControl = control,
            method = "svmRadial",
            allowParallel = T
          )
          ### 输出结果
          svm_res <- predictors(results)

          dir.create(file.path(od,'svm_rfe'),showWarnings = F,recursive = T)

          write_tsv(tibble(svm_res), file = file.path(od, "svm_rfe/svm_res.txt"))
          # 绘制结果
          pdf(file.path(od, "svm_rfe/svm_res.pdf"), height = 4, width = 4)
          plot(results, type = c("g", "o"))
          dev.off()

          return(results)
        }
      )
      res = map(y,~ f(.x))
      names(res) = y
      # tryCatch(f(),error = \(e){
      #   message('运行失败，请检查。')
      # })
    }
  )
)



if (F) {
  xgene <- read_delim("/data/Separated_Users/wangcy/project/P04_F220809002_osteoporosis/一开二次重启/data/Cellular iron metabolism.txt") %>% pull(1)

  load("/Pub/Users/wangyk/project/Poroject/P04_F220809002_osteoporosis/results/0.prepare_data/GSE56815_GPL96.RData", verbose = T)

  my_ml_obj <- ml$new(
    exprs = GSE56815_GPL96$data_exprs[xgene, ],
    pheno = GSE56815_GPL96$group %>% rownames_to_column("sample"),
    od = "/Pub/Users/wangyk/project/Poroject/F230908001_mr/test/"
  )

  my_ml_obj$help()

  a <- my_ml_obj$run("lasso")
  a <- my_ml_obj$run("svm")
  a <- my_ml_obj$run("rf")
  a
  class(a)
}


if (F) {
  suppressMessages(library(tidyverse))
  suppressMessages(library(magrittr))
  library(R6)

  as.integer(sqrt(22))

  ma <- R6::R6Class(
    classname = "ma",
    public = list(
      name = NULL,
      type = NULL,
      initialize = function(name, type, tui, ...) {
        self$name <- name
        self$type <- type
        private$tui <- tui
      },
      shangan = function(maan, ...) {
        self$maan <- maan
        paste0("gei ", self$type, " ", self$name, " shangge ", maan) %>% print()
        invisible(self)
      },
      chuanxie = function(...){
        print(paste('xuyao', as.integer(private$tui/2) ,'shuang xie.'))
        invisible(self)
      }
    ),
    private = list(
      tui = NA
    ),
    lock_objects = F
  )

  ma$set("public", "p", function(...) {
    paste(self$name, self$type)
  })

  zhenzhuxiaer <- ma$new(name = "zhenzhu", type = "xiaer", tui = 4)
  zhenzhuxiaer <- zhenzhuxiaer$shangan(maan = "huangjin")
  zhenzhuxiaer <- zhenzhuxiaer$chuanxie()
  zhenzhuxiaer$p()
}


# v2
if(F){
  ml <- R6::R6Class(
  classname = "ml",
  public = list(
    exprs = NULL,
    pheno = NULL,
    od = NULL,
    log = NULL,
    initialize = \(exprs = NULL, pheno = NULL, od = NULL, log = NULL){
      self$exprs <- exprs
      self$pheno <- pheno
      self$od <- od

      if (is.null(log)) {
        self$log <- logging$new("ml")
      }

      index <- list(self$exprs, self$pheno, self$od)

      if (any(map_lgl(index, is.null))) {
        features <- c("expr", "pheno", "od")[which(map_lgl(index, is.null))]
        self$log$warning("缺少{features}信息，需要在`new()`中重新定义")
      }
    },
    help = \(...){
      message("
      help文档：
        - 随机森林中，包含ntree, var_name，seed参数，第一个对应随机森林分析中的ntree数目，默认1000，第二个参数为生成结果的标识字段，默认为'randomforest'，seed 对应分析中的钟子号；
        - lasso中包含seed参数，默认1110;
        - svm-rfe默认使用五倍交叉验证；
        - 使用示例：my_ml_obj$run('lasso',seed = 123.1234)
      ")
    },
    run = function(method = NULL, ...) {
      if (is.null(method)) {
        self$log$error("请选择正确的参数，程序暂时不运行任何分析，参考{.arg method = 'lasso|rf|svm_rfe}'")
        self$help()
        return()
      }

      y <- match.arg(method, c("lasso", "rf", "svm_rfe"))

      f <- function(y) {
        ml_func <- switch(y,
          lasso = \(exprs = self$exprs, od = self$od, pheno = self$pheno, seed = 1110, ...){
            # lasso 建模
            library(glmnet)
            set.seed(seed)
            if (any(is.na(exprs[[2]]))) {
              self$log$warning("LASSO: 表达谱有NA，自动去除...")
              exprs <- na.omit(exprs)
            }

            x_tmp <- exprs %>%
              t() %>%
              as.data.frame() %>%
              rownames_to_column("sample")
            dat <- pheno %>% inner_join(x_tmp)

            x <- dat %>%
              dplyr::select(-any_of(colnames(pheno))) %>%
              as.matrix()
            y <- dat %>%
              dplyr::select(any_of(colnames(pheno))) %>%
              select(-sample) %>%
              pull(1)
            fit1 <<- cv.glmnet(x, y, family = "binomial", type.measure = "deviance", alpha = 1)
            # 绘制图片1
            dir.create(file.path(od), recursive = T, showWarnings = F)
            pdf(file = file.path(od, "lasso_cv.glmnet.fit.pdf"), width = 5, height = 4)
            plot(fit1)
            dev.off()
            # 查看模型各变量系数
            fit1_coef_lambda.min <- coef(fit1, s = fit1$lambda.min)
            # 查看系数非0变量
            fit1.min <- fit1_coef_lambda.min[which(fit1_coef_lambda.min != 0), ]
            # 矩阵转换，使上面结果更方便查看
            fit1.min_2 <- matrix(fit1.min, length(fit1.min), 1)
            rownames(fit1.min_2) <- names(fit1.min)
            colnames(fit1.min_2) <- c("coef")
            lasso_res <- fit1.min_2 %>%
              as.data.frame() %>%
              rownames_to_column("gene")
            lasso_res <- lasso_res[-1, ]
            # (Intercept) -98.7319020

            # 绘制图片2
            fit2 <<- glmnet(x, y, family = "binomial")

            pdf(file = file.path(od, "lasso_glmnet.fit.pdf"), width = 5, height = 4)
            plot(fit2, xvar = "lambda", label = TRUE)
            abline(v = log(fit1$lambda.min), col = "black", lty = 2)
            dev.off()

            lasso_fits <- list(x = x, y = y, fit = fit2, cv.fit = fit1)
            saveRDS(lasso_fits, file.path(od, "lasso_fits.rds"))

            # 输出结果
            write_tsv(lasso_res, file.path(od, "lasso_res.txt"))

            self$log$info("LASSO binomial 回归分析结束")
            lasso_fits$lasso_res <- lasso_res
            lasso_fits
          },
          rf = function(data = self$exprs, group = self$pheno, od = self$od, 
          ntree = 1000, var_name = "random_forest",seed = 1110, ...) {
            self$log$info("随机森林分析开始")
            data <- na.omit(data)

            data <- data %>%
              t() %>%
              as.data.frame() %>%
              rownames_to_column("sample") %>%
              inner_join(group)

            group_col <- colnames(group)[2]

            data[[group_col]] <- factor(data[[group_col]])
            data <- data[, -1]

            dir.create(od, recursive = T, showWarnings = F)

            colnames(data) <- colnames(data) %>%
              str_replace_all("-", "_")

            library(randomForest)
            library(caret)
            set.seed(seed)

            mtry_range <- 1:(ncol(data) / 2)

            # 定义训练控制参数，使用 5 折交叉验证
            train_control <- trainControl(method = "cv", number = 5)

            # 定义 mtry 的候选值
            tune_grid <- expand.grid(.mtry = mtry_range)

            self$log$info("5折交叉验证进行mtry超参数调优")
            rf_model <- train(
              y = data[[group_col]],
              x = data %>%
                select(-any_of(group_col)),
              data = data,
              method = "rf",
              trControl = train_control,
              tuneGrid = tune_grid, importance = T, ntree = ntree
            )

            # bestTune = rf_model$bestTune
            # self$log$info("mtry最优数值：{bestTune}")

            a <- ggplot(rf_model)
            a$data$type <- "normal"
            a$data$type[which.max(a$data[[1]])[1]] <- "best"

            p <- a + geom_point(aes(color = type), show.legend = F) +
              scale_color_manual(values = c("best" = "#ca3434")) +
              egg::theme_article()

            pdf(paste0(od, "/01.", var_name, "_mtry准确率率波动曲线.pdf"), 4.5, 4)
            print(p)
            dev.off()

            pdf(paste0(od, "/02.", var_name, "_gini_准确度.pdf"), 6, 7.4)
            varImpPlot(rf_model$finalModel, main = NULL)
            dev.off()


            pdf(paste0(od, "/03.", var_name, "_错误率波动曲线.pdf"), 5.6, 4.5)
            layout(matrix(c(1, 2), nrow = 1),
              width = c(4, 1)
            )
            par(mar = c(5, 4, 4, 0)) # No margin on the right side
            plot(rf_model$finalModel, main = NULL,lwd= 2,lty=1)
            par(mar = c(5, 0, 4, 2)) # No margin on the left side
            plot(c(0, 1), type = "n", axes = F, xlab = "", ylab = "")
            legend("top", colnames(rf_model$finalModel$err.rate), col = 1:4, cex = 0.8, fill = 1:4)
            dev.off()

            # rf <- randomForest(
            #     y = data[[group_col]], x = data %>%
            #         select(-any_of(group_col)), importance = TRUE, proximity = TRUE, ntree = ntree,
            #     mtry = rf_model$bestTune
            # )
            d <- as.data.frame(importance(rf_model$finalModel)) %>% arrange(desc(MeanDecreaseGini))
            saveRDS(rf_model, file = paste0(od, "/04.", var_name, "_随机森林分析结果.rds"))
            saveRDS(d, file = paste0(od, "/05.", var_name, "_随机森林分析结果_importance_df.rds"))


            self$log$info("随机森林分析结束")
            return(list("d" = d, "train_rf_obj" = rf_model, "best_mtry" = rf_model$bestTune))
          },
          svm_rfe = \(exprs = self$exprs, od = self$od, pheno = self$pheno, seed = 111213,...){
            # SVM
            set.seed(seed)
            features <- exprs %>%
              t() %>%
              as_tibble(rownames = "sample")

            labels <- pheno %>% as_tibble()

            exp_t <- inner_join(features, labels, by = "sample") %>%
              column_to_rownames("sample")

            exp_t[is.na(exp_t)] <- 0

            x <- as.matrix(exp_t[, 1:(ncol(exp_t) - 1)])
            # x <- exp_t
            y <- factor(exp_t[, ncol(exp_t)])

            library(caret)
            library(kernlab)

            self$log$info("SVM-RFE 10x cross validate...")

            control <- rfeControl(functions = caretFuncs, method = "cv", number = 5)
            results <<- rfe(as.matrix(x),
              y = y,
              sizes = seq(1, ncol(x)/2),
              rfeControl = control,
              method = "svmRadial",
              allowParallel = T
            )
            ### 输出结果
            svm_res <- predictors(results)

            self$log$info("SVM-RFE 结束")

            dir.create(od,recursive = F,showWarnings = F)

            write_tsv(tibble(svm_res), file = file.path(od, "svm_res.txt"))
            # 绘制结果
            pdf(file.path(od, "svm_res.pdf"), height = 4, width = 4)
            plot(results, type = c("g", "o"))
            dev.off()

            saveRDS(results,file.path(od, "svm_results.rds"))

            return(results)
          }
        )
        return(ml_func)
      }

      do.call(f(y), list(...))
      # tryCatch(f(),error = \(e){
      #   message('运行失败，请检查。')
      # })
    }
  )
)

}